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Abstract 

A Keller-Segel model describes macroscopic dynamics of bacterial colonies and 
biological cells as well as dynamics of a gas of self-gravitating Brownian parti- 
cles. Bacteria secret chemical which attracts other bacteria so that they move 
towards chemical gradient creating nonlocal attraction between bacteria. If bac- 
terial (or Brownian particle) density exceeds a critical value then the density 
collapses (blows up) in a finite time which corresponds to bacterial aggregation 
or gravitational collapse. Collapse in the Keller-Segel model has striking qualita- 
tive similarities with a nonlinear Schrodinger equation including critical collapse 
in two dimensions and supercritical collapse in three dimensions. A self-similar 
solution near blow up point is studied in the critical two-dimensional case and it 
has a form of a rescaled steady state solution which contains a critical number 
of bacteria. Time dependence of scaling of that solution has square root scaling 
law with logarithmic modification. 
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1. Introduction 

Formation of singularity in a finite time (blow up) is a quite general phe- 
nomenon in many nonlinear systems including self-focusing in nonlinear optics, 
plasmas, hydrodynamics, and collapse of Bose-Einstein condensate [H, 0|- Blow 
up is often accompanied by a dramatic contraction of the spatial extent of 
solution, which is called collapse Collapse typically occurs when there is 
self-attraction in a nonlinear system. Here a collapse in a system of bacteria or 
biological cells is considered. Below we refer to bacteria and cell as synonyms. 
Bacteria often communicate through chemotaxis, when bacteria move towards 
chemical gradient of substance called chemoattractant. Bacteria can also se- 
crete the same chemoattractant which creates nonlocal attraction between them 
through secretion, diffusion and detection of chemoattractant by other bacteria 
of the same type. Bacteria are self-propelled and without chemotactic clue cen- 
ter of mass of each bacteria typically experiences random walk. Motion of bac- 
terial colonies is thus affected by permanent competition of random- walk-based 
diffusion and chemotaxis-based attraction. Macroscopically averaged motion of 
bacteria can be described by a Keller-Segel model (see e.g. [3-17] and references 
therein) : 

dtp = DV^p-V[kpVc], (1) 
dtc = Dc\/^c + ap, (2) 

where p(r, t) is the density of at spatial point r and time t, c is the concentration 
of chemoattractant, D is the diffusion coefficient of bacteria due to random walk, 
the coefficient fc > characterize strength of chemotaxis, Dc is the diffusion 
coefficient of chemoattractant and a is the production rate of chemoattractant 
by bacteria. 

If density of bacteria is initially low then diffusion typically exceeds attrac- 
tion and density remains low. E.g., a typical time scale for evolution of the a 
low density bacteria Escherichia coli distribution in a petri dish is about a day 
d (see Fi gure 3 A in Ref. Q). If initial density is not very small so that attrac- 
tion wins over diffusion then bacteria experience aggregation (see Figure 3B in 
Ref. Q). A typical time scale of such aggregation in experiment on Escherichia 
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coli [3| is several minutes so it has an explosive character compare to evolution 
of bacteria outside of aggregation area. A goal of this Letter is to study how 
aggregation occurs from moderate initial bacterial density (below we give exact 
criterion for formation of aggregation) . Aggregation corresponds to the collapse 
of bacterial density in the Keller-Segel model ([T|), Aggregation is a first 
step to a formation of multicellular organisms and quite important in biological 
applications Q. Near singularity the Keller-Segel model is not applicable when 
typical distance between bacteria is about size of bacteria. In that regime mod- 
ification of the Keller-Segel equation was derived from microscopic stochastic 
dynamics of bacteria which prevents collapse due to excluded volume constraint 
(different bacteria cannot occupy the same volume) [l8l.[l9l.[20t. Here however 



the original Keller-Segel model without regularization is considered. 

Typically diffusion of a chemical is much faster than cell diffusion, i.e. 
p/Dc < 1. E.g., D/Dc ^ 1/40- 1/400 for the cellular slime mold Dictyostelium 
[iil, and D/Dc 1/30 for microgha cells and neutrophils [H, H^]. Thus Eq. 
^ evolves on a much smaller time scale than Eq. ([T]) and one can neglect time 
derivative in 1^ which gives a reduced Keller-Segel equation (RKS) 

ftp- VV-V[pVc], 
V^c = -p, 

where we assumed that D, Dc,a and k are constants and rewrote all variables 
in dimensionless form as i — ?► tot, r — )• tl^^D^^'^r, p — > {Dc/toCxk)p, c — > {D/k)c 
and to is a typical timescale of the dynamics of p in Eq. ([T|) . 

Eq. ^ also describes a dynamics of a gas of self-gravitating Brownian 



particles which has applications in astrophysics [12|, |17[ . In that case the second 
Eq. in ([3]) is a Poisson equation for the gravitational potential — c while p is the 
gas density (all units are dimensionless) . The first Eq. in (jS]) is a Smoluchowski 
equation for p in the gravitational potential — c. Below we refer to p and c as 
the density of bacteria and the concentration of chemoattractant, respectively, 
but all results of this Letter are equally true for a gravitational collapse of the 
gas of self-gravitating Brownian particles. 

The main result of this work is that solution of RKS ([3]) in a spatial dimension 
D = 2 near collapse time to in neighborhood of spatial location of collapse r = 
has the following cylindrically-symmetrical self-similar form 

1 8 

P 



L{tY (1 + 2/2)2^ 
c = -21n(l + y2), 



y 



L{ty 



(4) 



L{t) = 2e ~y/to^e 



^ ^ ^( {\n[-\n{to-t)]f 
\ ^-2\n{to-t) 



where 7 — 0.577216 ... is the Euler's constant and L{t) = L is the time- 
dependent spatial width of solution (below we often omit argument of L for 
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brevity). This result has striking similarities with critical collapse in a nonlin- 
ear Schrodinger equation (NLS) [24-33]. E.g., collapsing solution of NLS has a 
form of a rescaled ground state soliton solution of NLS while (|1]) has a form of a 
rescaled steady state solution of ([3]) . Total number of bacteria in a self-similar 
solution (U) iV = J pdr = 8w is conserved which is similar to conservation of the 
total number of particles in collapsing NLS solution. Scaling of collapsing NLS 
solution has the same leading order form (to — t)^^"^ as in Q but a correction for 
NLS has a well-known log-log form L(t) = (27r)i/2(to-t)i/V(ln [- In (<o - t)])^/^ 
which is different from an exponent-of-square-root-of-logarithm for @. The 
logarithmic corrections in both RKS and NLS result from nearly exact balance 
between linear and nonlinear terms (respectively balance between diffusion and 
attraction for RKS and balance between dispersion and nonlinear attraction in 
NLS). The log- log scaling for NLS is quite challenging for numerical verifica- 
tion (sij. E.g., it is difficult to distinguish numerically the log-log scaling from 
a power of log scaling. In contrast, the exponent-of-square-root-of-logarithm 
scaling Q is much faster and gives sizable corrections even for a moderate in- 
crease of an amplitude of p. Similar scaling to dH) was obtained in Refs. 0, [U 
using a formal matched asymptotic expansion of RKS. The order of error in 
(jH) is however quite different from both and 11 1 because of power 3 of the 
log-log factor. In Ref. [llj a power of the log-log term is 1 instead. Result 

of Ref. predicts L{t) cx y^t^e-^^'^i- In {to - t)](^/*^(- inCo-t))-^'" 
scaling which only partially agrees with ([4]). This difference is not a big sur- 
prise because the formal matched asymptotic expansion of Refs. [7^ and [ll| 
does not provide a good control of errors. Another competing scaling is L{t) oc 



y/tn — te 2 a/ 2 which was obtained in Ref. [ij] (see also Ref. 



■(to-*)] 



17[ for more discussion on that scaling). It remains a challenge for a future 



numerical simulations to verify which of these predictions of L(t) are correct. 
Previous simulations 13, 12 1 showed existence of corrections to the leading order 
scaling L{t) cx \/to^^ but were not able to resolve the logarithmic corrections. 

Principle difference of the analysis of RKS ([3]) and derivation of the scaling 
^ in this Letter from all previous works is the use of ideas from the NLS collapse 
analysis. It allows to find of an appropriate gauge transformation which brings 
a linearization operator about solution (|4]) to the self-adjoint form. 

In Section [2] we study general properties of collapse of RKS, prove collapse 
existence and explore multiple similarities with collapse in NLS. Section |3] fo- 
cuses on a self-similar solution of 2D RKS. We write the self-similar solution as 
a rescaled steady state solution in blow up variables. We choose a gauge trans- 
formation to a new dependent variable for perturbations about the self-similar 
solution to bring a linearization operator about the self-similar solution to a self- 
adjoint form. In section |4] we expand these perturbations into eigenfunctions of 
the self-adjoint linearization operator to derive a set of amplitude equations for 
coefficients of that expansion. Compatibility conditions to have an adiabatic 
form of that expansion result in a closed differential equation (|47p. In Section 
[5] we solve the equation (|47| to derive the scaling Q. In Appendix we analyze 
a spectrum of the linearization operator and approximate its eigenfunctions. 
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2. Collapse in the Reduced Keller-Segel Equation 

RKS ([3]) can be written in a form of a conservation law 

dtp = -v-r, 

where F is the flux of bacterial probability density given by 



In p 



(5) 



(6) 



and c(r) is determined by the fundamental solution E{r, r') of the Poisson equa- 
tion. Here we consider a two-dimensional bacterial motion with 



;(r) = - / E{r, r')p{r')dr', E{r, r') = -!- In |r - 
J 27r 



(7) 



Eq. dl]) allows to rewrite RKS ([3]) as a closed integro-differential equation for p. 
Nonlocality of interaction is due to diffusion of chemoattractant. 

Assuming decaying boundary conditions at infinity we obtain a conservation 
of a total number N of bacteria: 



N = J p{r)dr = const. 
One can also define a Lyapunov functional 

p(r)c(r) 



£ = 



p{r) In p{r) — p{r) — 



dr, 



and represent RKS in a gradient form 



.S£\ 6E 



5p 



5tp = V-LoV— , —=\np-c, 



6p 



where the Lyapunov functional £ is a non-increasing function of time 

r2 



d£ 
It 



-dr. 



P 



(8) 



(9) 



(10) 



(11) 



£ is conserved only for a steady state solutions with zero flux F = 0. 

Although RKS ([3]) is a gradient non-Hamiltonian system, it has many strik- 
ing similarities with NLS [32,]: 



iStV' + V^V'-HlV'PV'^O. 



(12) 



NLS can be written in a Hamiltonian form idt'4' = with the Hamiltonian 



H = J jVV'P — dr. NLS has huge number of applications, e.g. it describes 
dynamics of Bose-Eistein condensate, self-focusing of light and propagation of 
nonlinear dispersive waves. Solutions of both RKS and NLS exist globally in 
dimension one D — 1. Dimension Z3 = 2 is critical for both RKS and NLS with 
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collapsing solutions existing for D > 2 [1, 24, 25, 2^ 34, 3^ 32 1. Collapse in 
critical dimension D = 2 is strong for both RKS and NLS, i.e. a finite number 
of bacteria (particles) is trapped into collapsing spatial region. For supercritical 
case D > 2, e.g. for D = 3, collapse in both RKS and NLS is weak which means 
that collapse is so fast that particles (bacteria) cannot keep up with collapse 
rate and close to collapse time t to a vanishing number of bacteria (particles) 
are trapped into collapsing region. 

To prove existence of collapse in RKS one can use a positive-definite quan- 
tity A = J r^pdr, which determines a mean square width of bacterial density 



distribution [36|, |37[ . Vanishing of A guarantees existence of collapse because of 
conservation of N. Proof of collapse existence in NLS ([T^ is based on a virial 
identity dfB — 8H [25|, |26[ and vanishing of B, where B = J r^jV'prfr is the 



positive-definite quantity. B is the analog of A in RKS because in quantum 
mechanical interpretation of NLS is the probability density of particle number, 
i.e. the analog of p in RKS. However, RKS is the non-Hamiltonian system 
and the direct analogy with a virial theorem for B does not work. Instead one 
can calculate time derivative of A using Eqs. ([3]), ([7]), integration by parts and 
vanishing boundary conditions at infinity which gives for D = 2 : 

■>^«^**'.«-f!, ,13, 



At=AN- 



±/2r.(r-r') 



where we also used symmetrization over r and r'. One concludes from that 
At < Q for N > Nc = 8tt and A turns negative in a finite time which proves 
existence of collapse. Existence of the critical number of bacteria Nc — 8tt is 
another similarity with NLS, where the critical number of particles Nc.nls = 
/ IV'pdr ~ 11.70. . . Qualitative difference between RKS and NLS ^ is that 
RKS is the integro-differential equation while NLS ([T2|) is a partial differential 
equation (PDE). However, it was shown in Ref. 38| that the virial identity can 



be used for a proof of collapse existence in an integro-differential equation of 
NLS-type which describes Bose- Einstein condensate with nonlocal dipole-dipole 
interaction. 



3. Self-similar collapsing solution of the 2D reduced Keller-Segel equa- 
tion 

Below we consider critical case D = 2 only. Both RKS and NLS are scale 
invariant in 2D which means that if p{r,t), c{r,t) is a solution of RKS then 
L^^p{r/L,t/L'^), c{r/L,t/L^) is also a solution for any L{t) = L = const > 0. 
Similar property holds for NLS. RKS has a static solution 

_ 8 

PO - (1+^2)2' (^4) 

Co = -21n(l + r2), 

with the critical number of bacteria N{pq) — 87r, which is another striking 
similarity with the ground state soliton solution of NLS ([T2|) also containing the 
critical number of particles Nc.nls- 
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Near collapse time t — > to and collapse location r = the RKS solution of 
p{r, t) approaches to a cylindrically-symmetric self-similar solution in a form of 
a rescaled stationary solution (|12p with time-dependent scale L : 



1 / r 



(15) 



where L{t) ^ for t — > to- 

We introduce a new dependent variable m for cylindrically-symmetric solu- 
tions: 

m(r,t) = -!-/ p{r',t)dr', (16) 



2tt 



r'\<r 



which allows 0] to rewrite RKS as a closed equation for m : 

dtm — rdrr^^drfn + r^^mdrfn. (17) 

m{r,t) in (I17|) has a meaning of a number of bacteria inside a circle of radius 
r (up to a factor 27r). Boundary condition for m at r — >■ oo is simply related 
to the total number of bacteria: m|r.=oo = N/{2tt). In contrast to RKS, pT)) is 
PDE for m. This simplification is possible for cylindrically-symmetric solutions 
of RKS only. 

The steady state solution (ITi)) of RKS takes the following form for m : 

"^o = Y^' (18) 
and the self-similar solution ((T51) becomes 

_ Jy^ 

^sel f similar ^ , o ; 

1 + y2 

y = f (19) 

Boundary condition at infinity gives the critical number of bacteria: 

27rmsei f similar —> 8tt — const which is also the indication that bacterial 

I y >oo 

collapse is strong because the number of bacteria trapped into collapse is con- 
stant. 

Assuming a power law dependence L{t) ^ (to — t)^ in the self-similar so- 
lution ([T9|) one concludes that all terms in Eq. p7|) are of the same order 
provided (3 — 1/2, which is similar to NLS. Like NLS, the self-similar solution 
(fTO|) is not an exact solution of Eq. p7)) and it is necessary to consider the 
logarithmic correction to L{t) ^ (to — t)^/^: L = {to ~ t)^/^/(ln (to — t)), where 
/(In (to — t)) is the slow function compare with (to — t)^/'^. That slow function 
comes from nearly exact balance between linear and nonlinear terms of RKS 
(between diffusion and attraction) and allows to introduce a small parameter 

a = -L{t)dtL{t) > 0, (20) 
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which is a slow function of {to — ty^"^ compare with L. A balance between linear 
and nonlinear terms of RKS becomes better with decrease of a — ?> . 

Based on the analogy with the critical NLS we introduce in Eq. pT)) new 
independent "blowup" variables (39j : 



dt' 



m 



(21) 



which produces from Eq. (1171) an equation for a new unknown function tp{y, r) = 
m{r, t) : 

drip = ydyU^^dyip + y^^tpdytp - aydyip, (22) 

where a is given by (|20|) . Near collapse L(t) — > which means that r — > oo, i.e. 
collapse occurs at infinity for the new time r. The advantage of moving to blow 
up variables is that the function Lp has bounded derivatives and the collapse 
time to is eliminated from consideration (it is mapped to r — )■ oo). 

Because we assume that a is a slow function, it is natural to look at solutions 
of Eq. (|22p in an adiabatic approximation when one can neglect r-derivative 
in the left hand side (Ihs) of Eq. (P^ . Then, assuming that \a\ <C 1, one 
can expand a solution of (|22l) in powers of a starting from ([T^ in the power 
zero. However, term —aydyf grows with y and violates that expansion for large 
y. So the adiabatic approximation can only work locally for not very large y 
which is also familiar in analysis of collapse in NLS. This however does not 
create big problems because for t — >■ it is important to take care of what is 
happening in the neighborhood of collapse point y = where the term —aydyip 
is indeed a small correction. Behavior of the solution at large y does not affect 
the self-similar solution near a collapse time. 

It is convenient to present a general solution of Eq. (1^^ in the following 
form 

4;/^ , s^y2 y^ 
1 + + 1 

where v{y,T) is responsible for all corrections with respect to the self-similar 
solution (fT9|) . Substitution of (fSSj) into ([22l) gives the following Eq.: 

drv + LaV = F. (24) 

Here 

U = --^dyy^'dy - + [^y' -2a + . (25) 

is a linear operator corresponding to a linearization of (|22l) with respect to ([TO]) , 
and the term 

F - - ^° c-a^'li + ^f2^^ay-/i , ^Z^^ ,2^, yvv' . 

^y y^ + l ^2(j/2 + l)^ +(^2 + 1)2^ +y2 + l^ 

(26) 



^{y, t) = + e^^y -i^v{y, r), (23) 
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is responsible for all other terms including nonlinear terms in v, an inhomo- 
geneous and linear terms. Generally F cannot be zero because (|19p is not an 
exact solution of for nonzero a. Eqs. are equivalent to Eq. (|22l) 

because up to now we have not made any approximations. 

Advantage of the definition (^5]) is that the operator Ca = —-^dyy^dy + V{y) 
has a form of a radially symmetric Schrodinger operator in spatial dimension 
four {D ~ 4) with a potential 

It means that Ca is the self-adjoint operator with a scalar product 

oo 

(7^,0) = J i^[v)<i>{y)y^dy. (28) 



The potential — > oo for y — > cx) which ensures that Ca has only discrete 
spectrum. This allows to expand arbitrary u in a discrete set of eigenfunctions 

of Ca. 

V = Citpi + C2'02 + C3tp3 + . . . , (29) 

where tpj (y) are the eigenfunctions of Ca , 

Calpj = Xjipj (30) 

and \j arc the respective eigenvalues. 

Note that a use of the scalar product (|28p which corresponds to the radially 
symmetric Schrodinger operator in D = 4 is simply the auxiliary mathematical 
trick because the operator Ca is self-adjoint with that scalar product. However, 
all solutions obtained below correspond to RKS ([3]) with D = 2. 

4. Amplitude equations 

In this Section we derive a set of amplitude equations for ci , C2 , . . . from 
(|29|) which solve Eq. (|24)) . We solve these amplitudes equations exploiting the 
fact that at the leading order in a the solution of p2|) is given by (fT9|) which we 
used in definition of (|23p . We expand all expressions below in integer powers of 
small parameters a and keeping a lowest nontrivial order of a. 

We show in Appendix that one can approximate three lowest eigenfunctions 
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ipi, ip2 and ip3 of Ca as follows: 



'02 
03 



1 + r 



1 + r 



1 



ay 

2 



l+J/2 

+a2/2in(l + y2) 



^ln(l + .^)je--V4, 

2 



12 + 7r2(2 + 7-ln2) 



12 



12 



7r2- 12 



In- 



12 



3 - 7 + In 2 



24 



12 



(31) 



where 7 = 0.577216 ... is the Euler's constant. The discrete spectrum of Ca is 
given by (see Appendix and Ref. (ioj): 



Ai=a(-2 + -^ + 2(l + 7-ln2)— ^ 



A2 
A3 



a ^ a' 



o 



(Ini 



1^3 



+ 0(a2) 



0(a2), 



(32) 



a 2 



Ini 



(Ini 



1\2 



+ 0(«'), 



and for j > 3 the eigenvalues are given by Eq. (jA.19|) . 

We assume (based e.g. on numerical simulations in [1, [l3|) that a is the 
adiabatically slow function of t: Cr ^ a^. An expansion of all amplitude 
equations in integer powers of a and -j^^ leads to a conclusion that drO, ^ 

a20(j^). We introduce a normalized function Ur — dra/a? = 0(1) which is 

oc drd and depends on In ^ only. A subscript t indicates that comes from 
drd- This allows to write drO- as an expansion in inverse powers of In 



drCL = 



Ini 



" In- 



O 



(33) 



where the coefficients a'r^ and Ot^'' are 0(1) and do not depend on r in the 
adiabatic approximation. 

Assume that initially the expansion coefficients Ci, C2, C3,... in are 
0(1). We substitute the expansion ((29|) of u from (|23p into (|24|) and scalar 
multiply ([M]) onto V'lj 02 e.t.c. (with the scalar product (^5)) ) to obtain a set 
of amplitude equations for ci,C2,.... A series expansion of these amplitude 



10 



equations into small a using Eq. p3p results in 

drCi + a — 2aci + O ^]~~T^ ~ 0' 



5^C3 + 2ac3 + O ( ^ ) = 



« ^ _ n (34) 



drCi + 4ac4 + O 



a 



where the terms 2a(j — 2), j = 1, 2, 3, 4, . . . originates from eigenvalues for ij^j, 
respectively (see Eq. (lA.lQp ). the term a in the first equation comes from the 
scalar product of -01 with the second term in the right hand side (rhs) of (|26p 
and we used the orthogonality of the vectors V'l, "02, "^Sj • ■ • ■ One concludes from 
Eqs. (j34p that the coefficients C3, C4, . . . initially decay exponentially until they 
reach the quasi-steady (adiabatic) state with 

C3,C4,...-0('^). (35) 

In the first equation of ([34|) we assume that 

«4 + 0(j2j) (36) 

to avoid exponential growth of ci in r. 

We have now a freedom in C2 and we choose it to have 1; — >■ for any y as 
a -> 0. According to (l5Tt . '0i(2/)ly=o = "02(2/) I y=o = 8 so we have to set 

^ a ^ 

In that case ci-0i + 02')/' = 0(a) for y = 0(1), i.e. v from vanishes with 
vanishing of a. 

Eqs. (1351) .po p and p7p justify the adiabatic approximation which means that 
the coefficients ci, C2, C3, C4, . . . depend on r through a only and one can series 
expand them in inverse powers of In ^ : 

i + 4"iii + <if'4^ + o(^)+o(„), 

1 7(1) 1 , ^(2) 1 " 1 



C2 = -77 + 4 TT^+O +0(a), 

2 ' (lni)2 V(lni)V (38) 



7(1) 1 ^^(2)_ 1 -I ^ 



11 



where the expansion coefficients dp-* = 0(1) for any i, j and they do not depend 
on T in the adiabatic approximation. 
It follows from IMl) and (El) that 



drc,=0[drj^ ) =0 



J = 1,2,3,, 



(39) 



Similar to Eqs. (1341) one can write amplitude Eqs. from the scalar multiplica- 
tion of dMl) by V'l, '4'2 and respectively. Using (IMl)-(Iini),(Ea)-(ISIl),(ISSl),(I3Hl), 
and series expanding all expressions in inverse powers of In - give the following 
equations: 



drCi 



a 
h^i 



~(o) 



(Ini 



l\2 



2d 



(1) 



2d 



(2) 



2d 



(1) 



O 



= 0, (40) 



OtC2 + f 

in T 



4 - 



(Ini 



1^2 



-l-f^ + 2dW 



+4°)(4'^ -24'V7 + ln2) 



5rC3 + r 

in - 



2d 



(1) 



(Ini 



l\2 



= 0. 



0, (41) 
(42) 



Here terms drCi, drC2 and drC3 can be omitted because they are O (^j^^^js^ 

according to Eq. ([M)) . Eqs. (P(I)) - (H^ are nothing more than compatibility 
conditions which ensure that expansions ([55)1 and ([55]) are correct so that a is 
indeed the adiabatically slow variable. 

It follows immediately from Eq. (^^ in the order that 



= 



(43) 



and from Eq. (pij) in the order that 



(44) 



and then from Eq. (HO]) in the order -j-S^ we obtain 



d« = -1. 

1 2 

Using Eqs. (gl]) and (|131)-(|l5]) we obtain in the order 



that 



(45) 



ai^^ = -2- 27 + 2 In 2. 



(46) 
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Eqs. and result in a closed differential equation for a: 



ar_ 2_ M / 1 \ 

lni + (lni)2+ Ulni)V' (47) 
M = -2-27 + 21n2. 

It follows from Eqs. ([55]) . P5)) and P5|) that the expansion coefficients are 
given by 

so that the coefBcients C3, C4, . . ., which correspond to positive eigenvalues, are 
of a lower order compare with ci and 02- One concludes that the self-similar 
solution (|19p is stable with respect to radially-symmetric perturbations and the 
leading order corrections to it are determined by a linear combination w ~ ciV'i + 
C2^2j where ci and C2 are given by (j48p . Typical dynamics of perturbations 
around the self-similar solution is shown in Figure 4 of Ref. [loj . 



5. Blow-up rate of self-similar solution 



In this section we solve Eq. (|T7)) to derive a blow-up rate (U). Integration of 
Eq. gZl) gives 



if mi 

a \ a 



M 

1 

2 



O 



1 



= 2t, 



(49) 



which is a version of the Lambert W-function and can be solved assuming r ^ 1 
by iterations for a as follows: 



1 lnln2T #-1 

In - =ln2T-lnln2r+— — + O 

a In 2t In 2t 



(lnln2r) 
(ln2T)2 



(50) 



Now we use Eqs. pO)) and ([2T|) to write a — —LdtL = —{drL)/L and express 
L through a as 



hiL = y a(T')dT', 



(51) 



where the lower limit of integration is irrelevant for r ^ 1. Using ([50)) and (|5ip 
we obtain that 



\nL = 



(ln2T)2 (ln2r)lnln2T M 



^ln2r + 0((lnln2T)3) 



(52) 
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We solve Eq. ^ ion t a.s to - t = L{tf{T')dT', where Lir) is given by 
52l) which results for r ^ 1 in 



-y + glng - y g + g + o({\nqf'^ 



to-t= (1/2) [q-' + 0{q-^ Inq)] exp 
q = ln2r, 



(53) 



where to is the collapse time. Solving Eq. (|53p by iterations for r ^ 1 we find 
an explicit expression for the rescaled time 



= 2- V-21n {to - i)eV-2in {to-t)+2+y 



1 + 



(ln[-ln(^o-t)])- 
v/-21n(to-i) 



(54) 



and using \/L{tY — dtT we obtain the scaling of the self-similar solution ([T^ : 



L(i) = 2e - Vto - te 



. / ln(tO-*)" 

- ^P+V 2 



1 + 



(ln[-ln(to-t)])^ 
^-21n(io-<) 



(55) 



A similar scaling was obtained in Refs. 0, [ll| using a formal matched asymp- 
totic expansion of RKS. Using Eq. ()55p and returning from the variable m of 
Eqs. (HI]) and (HT]) to the variables p and c of RKS © we obtain Eqs. (g]) which 
is the main result of this work. 



Appendix A. Eigenvalues and eigenfunctions of the linearization op- 
erator Ca 

The goal of this Appendix is to calculate the eigenvectors ipj and the eigen- 
values tpj of the operator £a It is impossible to write an explicit expressions 
for them through standard mathematical functions so we have to use approxi- 
mations. 

The operator Ca ([25l) is reduced to 

Co = -^dyy'dy ~ ^ - 2a, (A.l) 
yd (l + y2)^ 

for y <^ a~^^^ which is the Schrodinger operator for shallow potential well in 
D = 4,. Co corresponds to the boundary of disappearance of discrete spectrum 
(if we replace the potential in (|A.1|) by a potential — (^i^iyi , S > then discrete 

spectrum is absent). Eq. Cotp — Xip has only one (up to an arbitrary constant 
factor) bounded solution 

^0 = (A.2) 

1 + r 

which is a soft mode of that operator: Coipo — ^aipo because a <C 1. That 
soft mode corresponds to the breaking of scale invariance of the steady state 
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solution ([T5)) in Eq. (|17p because of a nontrivial time dependence of L{t) 



For the opposite limit y ^ a~^/^ the operator Ca (fSSj) is reduced to the 
the Schrodinger operator for the spherically symmetric harmonic potential in 
D = 4 : 

Ch = -^dyy'dy + ^y^~2a. (A.3) 

Solution of Ch'tp = AV' which decays for j/ — >■ cx) is given by 

V;,=ae-'^^^/4[/(-A,2,^), (A.4) 

where U{a, b, x) is the confluent hypergeometric function of the second kind j4l|. 
Here we use the factor a in rhs of (|A.4|) to match asymptotic of eigenfunctions 
for large y as will be seen below. An asymptotic of i/'/i for y S> a"-^/^ is given 



by m 



ijjh OC e-«yV4aA/2a+lyA/a 



l + O 



(A.5) 



Matching of perturbations about the asymptotics (jA.2p and (|A.4p allows to 
find the eigenvalues of Ca (US]) ^40|. In this work for calculation of the scalar 
products in Section |4] we need not only eigenvalues but eigenfunctions also. So 
instead of matching the asymptotics we use a variational approximation for the 
eigenvalues and the eigenvectors of Ca- Idea is to satisfy simultaneously the 
asymptotics (IA.2P and (|A.5P for optimal choice of trial eigenfunctions. For the 
first eigenvector -0i of Ca we choose 

V'l = T^e-"^'/4 (A.6) 

and approximate the eigenvalue Ai as 

A. = (A.7) 

Calculation of (jA.7[) assuming a ^ 1 gives 

Ai = « (-2 + ^ + 2(1 + 7 - 1-2)^) + O + 0{a^). (A.8) 

Here we have two small parameters a and In ^ . We keep a minimal number 
of terms which are necessary for derivation of L[t) dependence in Sections 0] and 
[5l One concludes from comparison of (lA.sp with (jA.6p that they have the same 
asymptotic tpi cx ?/~^e~°^ for y ^ a~^/^. One can now see a motivation of 
our choice of the extra factor a in (jA.4p because tpi in (IA.6P with such choice 
matches the asymptotics between y <^ a~^^^ and y 3> a^^^^ not only in power 
of y but also in power of a. 
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From an analogy with the spectrum of the spherically symmetric harmonic 
oscillator like (|A.3|) we expect that the eigenvalues of La ((25|) are separated 

at the leading order by 2a from each other. So that A2 — O (kf^) from 

(|A.5[) one concludes that '!/'2 « ae^"-^ 1^ for y ^ a^^l"^ . In addition, function 
U(a,b,x) has terms ocl/x and oc Inx for small x which all together motivates 
the following choice of a trial function for the second eigenfunction Tp2 '■ 

^2 = (1 + 4'^«y' + cp'ay^ In (1 + y')) e^^^/^, (A.9) 



y 

(2) (2) 

where the trial constants and C; should be chosen from an orthogonality 
condition 



(^Ai>2)=0 (A.IO) 
and from a minimization of the second eigenvalue 

A2 = mm — -— . (A.ll) 

cf\cf> {^2,^2) 

Our choice of the term In (1 + y'^) instead of In y^ in (|A.9I) was motivated by the 
necessity to avoid a logarithmic singularity for all real y. In addition, a similar 
term occurs if one approximates the eigenfunctions of Ca (|25p by a perturbation 
expansion near the soft mode (|A.2[) 40] . 

Eqs. (|A.9p and (jA.lOp result in the following expressions 



c. 



(2) _ 1 
2 ~ 2 ' 

h ^ 2 



(A.12) 



and 



A, . „ (-L + 2(2 + - 1„2,^) + O (^) + 0«.^). (A.13) 

The coefficients Cj^' and cp-* in (jA.9[) are chosen from a necessity to satisfy the 

orthogonality condition (jA.10[) up to the order o(^j^^. This choice is not 

unique because one can also use the minimization condition (lA.lip . However, 
use of that condition would result in change of values of A2 in (IA.13P in terms 



of O y j which is beyond an order we need to calculate Thus our 
motivation for the specific choice (|A.9p is to have the most compact expressions 
for Cj^-* and cp^ . It was checked that other possible choices of c^^ and cf^ do 
not change both A2 in (lA.lSp and final answer ^ because they correspond to 

higher order corrections O (^-fj^frjj^- 

From A3 ~ 2a and (|A.5P one concludes that -03 oc a^y^e""^ for y ^ a'^^"^. 
Together with (|A.9p that suggests the following choice for the third eigenfunction 

^3 = (1 + c^^^ay' + c!i^ay^ In (1 + y^) + c['\^y^) e-^^/^ (A.14) 
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where the trial constants c^^ , cf"^ and cf should be chosen from two orthogo- 
nality conditions 



.(3) 



and from a minimization of the third eigenvalue 



Eqs. (|A.14p and (|A.15P result in the following expressions 

TT^lni -12 + 7r2(2 + 7-ln2) 



(A.15) 



(A.16) 



.(3) 
-2 

„(3) 



7r2- 12 
12 



7r2 - 12 



TT^ - 12' 



(A.17) 



In- 



and 



2 

h^ 



3 - 7 + In 2 



O 



24 



7r2 - 12 



(A.18) 



Similar to ip2, the coefficients , cj and cjj in (jA.14l) are chosen from a 
necessity to satisfy the orthogonality conditions (jA.151) up to the order O ( i~t) • 
This choice is not unique because one can also use the minimization condition 
(|A.16p . However, use of that condition would result in change of values of A3 

in (IA.13[) in terms of O (^-^j^fr^^- But terms of that order are not needed to 

calculate (|4]) in Sections U and [5] 

Eqs. (US]), (TOl) . (IXT2]) . (1X131) . (IXTl) . (|XT7l) . (TOSj) resuh in Eqs. 

(PT|) and The eigenvalues (|5T|) are exactly the same as can be obtained 

from a matched asymptotic expansion [40| which ensures accuracy of amplitude 
equations in Sections 0] and [SJ As we mentioned above, choice of a trial functions 
is not unique. However we checked that different choices do not change results 
of Sections 2] and [5] because they correspond to the higher order corrections. We 
also confirmed accuracy of the approximations of the eigenfunctions pip and 
the eigenvalues p2l) by a direct numerical evaluation of the spectrum of £a (|25p 
for different values of a. 

Other eigenvalues have the following form (40| : 



+ 0(a2). 



(A.19) 



Corresponding eigenvectors can be obtained by an extension of the variational 
procedure described above through the inclusion of extra terms with the asymp- 
totic Ipj (X a^-l2/20-2)g-ayV4 foj. y ^ ^-1/2^ 
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